home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_13_06 / dwyer / randport.c < prev   
Text File  |  1995-02-10  |  6KB  |  226 lines
  1. /* Listing 2
  2.     rand_por[t].c
  3.     a portable and reasonably fast random number generator
  4.         multiplicative random number generator
  5.     see
  6.        L'Ecuyer - Comm. of the ACM, Oct. 1990, vol. 33.
  7.        Numerical Recipes in C, 2nd edition, pp.  278-86
  8.        L'Ecuyer and Cote, ACM Transactions on Mathematical Software,
  9.          March 1991
  10.        Russian peasant algorithm -- Knuth, vol. II, pp.  442-43
  11.     Copyright (c) 1994, 1995 by Gerald P. Dwyer, Jr.     */
  12.  
  13. #include    <time.h>
  14. #include    <stdlib.h>
  15. #include    <limits.h>
  16. #include    <assert.h>
  17.  
  18. #define    TESTING
  19.  
  20. #define    TRUE (-1)
  21. #define    FALSE 0
  22.  
  23. long init_rand_port(long seed) ;
  24. long get_init_rand_port(void) ;
  25. long genr_rand_port(long init_rand) ;
  26. long rand_port(void) ;
  27. double rand_rect_port(void) ;
  28. long skip_ahead(long a, long init_rand, long modulus, long skip) ;
  29. long mult_mod(long a, long x, long modulus) ;
  30.  
  31. #define  MOD    2147483647L     /* modulus for generator */
  32. #define  MULT        41358L     /* multiplier            */
  33.                           /* modulus = mult*quotient + remainder */
  34. #define  Q           51924L     /* int(modulus / multiplier) */
  35. #define  R           10855L     /* remainder                 */
  36. #define  MAX_VALUE    (MOD-1)
  37.  
  38. #define  EXP_VAL 1285562981L    /* value for 10,000th draw */
  39.  
  40. #define  IMPOSSIBLE_RAND (-1)
  41. #define  STARTUP_RANDS   16     /* throw away this number of
  42.                           initial random numbers */
  43.  
  44. static long rand_num = IMPOSSIBLE_RAND ;
  45.  
  46. /* initialize random number generator with seed */
  47. long init_rand_port(long seed)
  48. {
  49.     extern long rand_num ;
  50.     int i ;
  51.  
  52.     if (seed < 1 || seed > MAX_VALUE)        /* if seed out of range */
  53.         seed = get_init_rand_port() ;            /* get seed */
  54.  
  55.     rand_num = seed ;
  56.     for (i = 0; i < STARTUP_RANDS; i++)    /* and throw away */
  57.         rand_num = genr_rand_port(rand_num) ;    /* some initial ones */
  58.  
  59.     return seed ;
  60. }
  61.  
  62.  
  63. /* get a long initial seed for generator
  64.     assumes that rand returns a short integer */
  65. long get_init_rand_port(void)
  66. {
  67.     long seed ;
  68.  
  69.     srand((unsigned int)time(NULL)) ;    /* initialize system generator */
  70.     do {
  71.         seed  = ((long)rand())*rand() ;
  72.         seed += ((long)rand())*rand() ;
  73.     } while (seed > MAX_VALUE) ;
  74.  
  75.     assert (seed > 0) ;
  76.  
  77.     return seed ;
  78. }
  79.  
  80.  
  81. /* generate the next value in sequence from generator
  82.     uses approximate factoring
  83.     residue = (a * x) mod modulus
  84.             = a*x - [(a*x)/modulus]*modulus
  85.     where
  86.         [(a*x)/modulus] = integer less than or equal to (a*x)/modulus
  87.     approximate factoring avoids overflow associated with a*x and
  88.         uses equivalence of above with
  89.     residue = a * (x - q * k) - r * k + (k-k1) * modulus
  90.     where
  91.         modulus = a * q + r
  92.         q  = [modulus/a]
  93.         k  = [x/q]        (= [ax/aq])
  94.         k1 = [a*x/modulus]
  95.     assumes
  96.         a, m > 0
  97.         0 < init_rand < modulus
  98.         a * a <= modulus
  99.         [a*x/a*q]-[a*x/modulus] <= 1
  100.             (for only one addition of modulus below) */
  101. long genr_rand_port(long init_rand)
  102. {
  103.     long k, residue ;
  104.  
  105.     k = init_rand / Q ;
  106.     residue = MULT * (init_rand - Q * k) - R * k ;
  107.     if (residue < 0)
  108.         residue += MOD ;
  109.  
  110.     assert(residue >= 1 && residue <= MAX_VALUE) ;
  111.  
  112.     return residue ;
  113. }
  114.  
  115.  
  116. /* get a random number */
  117. long rand_port(void)
  118. {
  119.     extern long rand_num ;
  120.  
  121.                 /* if not initialized, do it now */
  122.     if (rand_num == IMPOSSIBLE_RAND) {
  123.         rand_num = 1 ;
  124.         init_rand_port(rand_num) ;
  125.     }
  126.  
  127.     rand_num = genr_rand_port(rand_num) ;
  128.  
  129.     return rand_num ;
  130. }
  131.  
  132.  
  133. /* generates a value on (0,1) with mean of .5
  134.     range of values is [1/(MAX_VALUE+1), MAX_VALUE/(MAX_VALUE+1)]
  135.     to get [0,1], use (double)(rand_port()-1)/(double)(MAX_VALUE-1) */
  136. double rand_rect_port(void)
  137. {
  138.     return (double)rand_port()/(double)(MAX_VALUE+1) ;
  139. }
  140.  
  141.  
  142. /* skip ahead in recursion
  143.     residue = (a^skip * init) mod modulus
  144.     Use Russian peasant algorithm          */
  145. long skip_ahead(long a, long init_rand, long modulus, long skip)
  146. {
  147.     long residue = 1 ;
  148.  
  149.     if (init_rand < 1 || init_rand > modulus-1 || skip < 0)
  150.         return -1 ;
  151.  
  152.     while (skip > 0) {
  153.         if (skip % 2)
  154.             residue = mult_mod(a, residue, modulus) ;
  155.         a = mult_mod(a, a, modulus) ;
  156.         skip >>= 1 ;
  157.     }
  158.     residue = mult_mod(residue, init_rand, modulus) ;
  159.  
  160.     assert(residue >= 1 && residue <= modulus-1) ;
  161.  
  162.     return residue ;
  163. }
  164.  
  165.  
  166. /* calculate residue = (a * x) mod modulus for arbitrary a and x
  167.         without overflow
  168.     assume 0 < a < modulus and 0 < x < modulus
  169.     use Russian peasant algorithm followed by approximate factoring */
  170. long mult_mod(long a, long x, long modulus)
  171. {
  172.     long q, r, k, residue ;
  173.  
  174.     residue = -modulus ;        /* to avoid overflow on addition */
  175.  
  176.     while (a > SHRT_MAX) {    /* use Russian Peasant to reduce a */
  177.         if (a % 2) {
  178.             residue += x ;
  179.             if (residue > 0)
  180.                 residue -= modulus ;
  181.         }
  182.         x += (x - modulus) ;
  183.         if (x < 0)
  184.             x += modulus ;
  185.         a >>= 1 ;
  186.     }
  187.                     /* now apply approximate factoring to a
  188.                         and compute (a * x) mod modulus */
  189.     q = modulus / a ;
  190.     r = modulus - a * q ;
  191.     k = x / q ;
  192.     x = a * (x - q * k) - r * k ;
  193.     while (x < 0)
  194.         x += modulus ;
  195.                     /* add result to residue and take mod */
  196.     residue += x ;
  197.     if (residue < 0)    /* undo initial subtraction if necessary */
  198.         residue += modulus ;
  199.  
  200.     assert(residue >= 1 && residue <= modulus-1) ;
  201.  
  202.     return residue ;
  203. }
  204.  
  205.  
  206. #if    defined(TESTING)
  207. /* Test the generator */
  208. #include    <stdio.h>
  209. void main(void)
  210. {
  211.     long seed ;
  212.     int i ;
  213.  
  214.     seed = init_rand_port(1) ;
  215.     printf("Seed for random number generator is %ld\n", seed) ;
  216.     i = STARTUP_RANDS ;   /* threw away STARTUP_RANDS */
  217.     do {
  218.         rand_port() ;
  219.         i++ ;
  220.     } while (i < 9999) ;
  221.  
  222.     printf("On draw 10000, random number should be %ld\n", EXP_VAL) ;
  223.     printf("On draw %d, random number is        %ld\n", i+1, rand_port()) ;
  224. }
  225. #endif    /* TESTING */
  226.